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ABSTRACT 

The  primary  objective  of  this  study  was  to  build  a 
mathematical  model  to  predict  the  probability  of  a  target 
moving  according  to  a  two-dimensional  random  tour  model 
avoiding  detection  (i.e.,  surviving)  to  some  specified 
time,  t. 


This  model  assumes  that  there  is  a  stationary  searcher 
having  a  "cookie-cutter"  sensor  located  in  the  center  of  the 
search  area. 
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A  Monte-Carlo  simulation^ program  was  used  to  generate 
the  non-detection  probabilities.  The  output  of  this  program 


was  used  to  construct  the  required  mathematical  model. 
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The  model  predicts,  and  simulation  supports,  that  as  the 
mean  segment  length  of  the  random  tour  becomes  small  with 
respect  to  the  square  root  of  the  area  size,  the  probability 
of  non-detection  approaches  that  previously  obtained  for  a 
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I.  RANDOM  TOUR  MODEL 


A.  INTRODUCTION 

The  main  objective  of  this  thesis  was  to  construct  and 
test  an  experimental  mathematical  model  to  predict  the 
probability  that  a  target  moving  according  to  a  two- 
dimensional  random  tour  will  avoid  detection  to  time  t  by  a 
fixed  sensor. 

B.  DESCRIPTION  OF  RANDOM  TOUR  MODEL 

1 .  The  Searcher  Location 
The  searcher  is  assumed  to  be  located  in  the  center 

of  a  square  search  region  of  area  A.  This  location  is  held 
fixed  during  the  search  period.  The  searcher  has  a  detec¬ 
tion  capability  over  a  disk  of  radius  R.  (See  Figure  1.1). 

The  detection  probability  for  a  target  inside  this 
disk  is  1  and  it  is  0  outside.  The  searcher  thus  has  a 
"cookie-cutter"  sensor  with  detection  range  R.  [Ref.  1] 

2.  The  Target  Starting  Position 

The  target's  starting  position  is  uniformly 
distributed  over  the  square  search  region  A. 

3.  Motion  of  the  Target 

The  target  moves  randomly  over  the  area  A  according 


to  a  random  tour  which  reflects  off  the  area  boundaries 


The  target  track  is  a  connected  sequence  of  line 
segments.  The  direction,  or  target  course  e,  for  each 
straight  segment  is  selected  from  an  independent  uniform 
distribution  between  0  and  2*  radians. 

The  length  of  time,  T,  the  target  spends  on  each  leg 
(assuming  no  reflection  off  the  area  boundaries)  is  selected 
from  an  independent  exponential  distribution  with  mean  1/x. 
The  term  x  is  the  rate  of  course  change  (again  ignoring 
reflections) . 

4 .  Detection 

Detection  occurs  the  first  time  the  target  enters 
the  searcher's  detection  disk;  that  is,  when  the  distance 
between  the  target  and  the  center  of  A  is  less  than  or  equal 
to  R. 


C.  NECESSITY  OF  SIMULATION 

An  analytic  expression  for  the  probability  density  of 
the  target's  position  after  a  random  tour  of  time  length  t 
was  derived  in  [Ref.  2].  Given  the  target's  initial 
position  at  the  origin  of  a  two  dimensional  coordinate 
system,  this  expression  is: 


g(x,y,t)  «  [e-xt/2*  (Vt)  2]  U(r-l)  +  [xt/  l[l-r2]exp(xt  jl-r2) } 

(1.1) 


where 


V  =  Target  speed  (nautical  miles  per  hour) 


*  =  Rate  of  course  change  (1/hour) 

t  *  Time  (hours) 


6  =  Dirac  «-function 

x,y  =  Components  of  the  target's  new  position. 

Expression  (1.1)  does  nr*:  account  for  boundary  effects 
and  it  considers  the  initiax  position  of  the  target  to  be 
the  origin.  Adding  the  effects  of  boundary  reflection  and 
assuming  the  initial  starting  position  to  be  uniformly 
distributed  over  A  significantly  complicates  the  calculation 
of  g  (x,y,t) . 

In  addition,  it  was  the  purpose  of  this  work  to  find  the 
probability  of  non-detection  to  time  t  (PND(t)),  not  the 
probability  density  function  for  the  target.  Thus,  it  was 
necessary  to  use  simulation  to  attack  this  problem. 

D.  SIMULATION  MODEL  OF  RANDOM  TOUR 

A  Monte-Carlo  simulation  computer  model  (called  Random 
Tour  Simulation  or  RATSIM)  was  used  to  estimate  PND(t)  for 
the  random  tour  model.  This  program  was  written  in  FORTRAN 
and  designed  to  run  on  the  IBM  3033  at  the  Naval 
Postgraduate  School.  It  uses  the  International  Mathematical 
and  Statistical  Library  (IMSL)  packages  GGUBS  to  generate 


uniform  random  variables  and  GGEXN  to  generate  exponential 


random  variables. 


1.  Inputs 


•  Radius  of  detection  disk  R,  in  nautical  miles  (nm) . 

•  Area  size  A,  in  square  nautical  miles  (nm  ). 

•  Target  speed  V,  in  nautical  mil*s  per  hour  (nm/hr). 

•  Rate  of  course  change  X,  in  1/hour  (hr  *) 

•  Number  of  replications  (REP). 

•  Detection  period  (TMAX) ,  in  hours  (hr) . 

•  Time  increment  AT,  in  minutes. 


2.  Functioning  of  the  Program 


(i)  At  the  start  of  each  replication,  the  initial 
starting  position  of  the  target  is  drawn  from  a 
uniform  distribution  over  the  area  A.  The  course  8 
is  drawn  from  a  uniform  distribution  on  (0,  2w). 

(ii)  The  course  is  changed  after  a  random  time  leg  T 
drawn  from  an  exponential  distribution  with  mean 
1/x. 

(iii)  After  each  time  increment  At,  the  new  position  of 
the  target  is  calculated  from: 


=  X  +  V  •  At  •  sin  0 


new  ~  old 


=  Y_i j  +  V  •  At  •  cos  e 


new  ~  old 


where 


Xnew'  *new  =  coordinates  of  the  new  position  at 


the  end  of  At. 


Xold'  *old  *  coordinates  of  the  old  position  at 
the  beginning  of  At. 


y+  * 
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Also,  the  distance  D  between  the  new  position  of  the 
target  and  the  center  of  the  searcher  disk  is 
calculated  from: 

d2  -  <x„e„  *  xser>2  +  <Vnew  '  xser>2  . 
where 


X  ,  Y  =  coordinates  of  the  center  of  the 
r  ser  searcher's  disk. 


(iv)  The  replication  terminates  if:  D  <  R  or  if  the 
detection  period  (TMAX)  is  over.  Then  a  new 
replication  begins.  The  process  continues  until  the 
specified  number  of  replications  is  reached. 

(v)  Two  counters  are  used,  one  to  determine  the  current 
time  t,  and  the  other  to  count  the  number  of 
replications  in  which  detection  occurs. 


3.  Design  of  the  Experiment 

Different  time  increments  At,  varying  from  1  minute 
up  to  10  minutes,  were  tested  with  RATSIM  and  3  minutes  was 
accepted  as  a  reasonable  compromise.  For  smaller  At,  the 
execution  time  of  the  program  increased  unacceptably.  For 
larger  values,  it  was  possible  for  the  simulated  path  to 
jump  across  a  significant  portion  of  the  detection  disk 
without  achieving  detection,  even  though  the  line  segment 
connecting  two  successive  discrete  positions  of  the  target 
was  partly  on  the  disk  [Ref.  3}.  This  will  reduce  the 
detection  rate,  especially  for  large  V.  However,  as 
illustrated  in  Figure  1.2,  PND(t)  can  be  relatively 
insensitive  to  At  less  than  10  minutes  when  the  problem 


Various 


parameters  are  appropriate  for  antisubmarine  warfare  (ASW) 
search. 

It  was  decided  to  conduct  2400  replications  for 
each  RATSIM  experiment.  This  resulted  in  the  standard 
deviation  of  the  simulated  PND (t)  being  no  greater  than 
0.25/2400  ■  0.0102. 

Also  the  maximum  time  allowed  for  detection  (TNAX) 
was  set  at  100  hours.  This  was  selected  to  allow  PND(TMAX) 
to  be  near  0  for  all  tested  values  of  problem  parameters. 

4.  Boundary  Effects 

When  the  target  encounters  a  boundary,  a  reflection 
is  made  to  keep  the  target  inside  the  search  area  A.  The 
target  position  after  reflection  is  determined  as  follows: 

In  Y-Direction: 

If  Y  <  0  then  Y  becomes  (-Y)  * 

If  Y  >  L  then  Y  becomes  (2L  -  Y)  . 

where 

L  *  length  of  a  side  of  the  square  search  area  A 
1  *e. ,  L  ■  A  . 

In  X-Direction: 

The  target  reflects  in  the  X-direction  in  a  similar 

manner. 

The  target  course  e  changes  after  reflection  as 


follows : 


At  Y  *  0  or  Y  =  L:  e  bcomes  (2*  -  e)  / 

At  X  =  0  or  X  =  L:  e  becomes  (*  -  0)  • 

Thus,  "the  angle  of  incidence  equals  the  angle  of 
reflection."  The  reflection  process  is  illustrated  in 
Figure  1.3. 

5.  Output 

At  each  time  t,  the  primary  simulation  output  is  the 
^T  -  N0 

ratio  — -  where 

T 

Ng  =  number  of  replications  giving  a  detection  by  time  t, 
and 

Nt  =  total  number  of  replications  used  in  Monte-Carlo 
simulation. 

This  ratio  is  the  simulated  probability  of  non¬ 
detection  by  time  t,  PND(t). 


II.  RELATIONSHIP  BETWEEN  THE  RANDOM  TOUR  AND 
DIFFUSION  MODELS 


A.  DESCRIPTION  OF  DIFFUSION  MODEL 

In  the  diffusion  model  considered  here,  the  target  moves 
randomly  over  a  square  search  area  A  according  to  Brownian 
motion  with  a  diffusion  constant  D  (units:  area/time). 


Perfect  reflection  occurs  at  the  area  boundaries. 


The  target  starting  position  is  uniformly  distributed 
over  A.  For  any  time  interval  of  length  At  which  does  not 
contain  a  boundary  reflection,  the  components  of  the 
target's  position  along  the  X  and  Y  axes  suffer  increments 
which  are  each  distributed  independently  and  normally  with 
mean  0  and  variance  D  At. 

A  searcher  having  a  "cookie-cutter"  sensor  with 
detection  range  R  is  located  at  the  center  of  the  search 


region. 


Detection  occurs  whenever  the  range  between  the  searcher 
and  the  target  becomes  R  or  less. 


B.  RELATIONSHIP  WITH  RANDOM  TOUR  MODEL 

In  [Ref.  4]  it  is  shown  that  as  the  rate  of  course 

change  x  for  an  unconstrained  random  tour  gets  larger  such 
2 

that  V  /x  ■  constant,  then  the  random  tour  can  be 
approximated  by  a  diffusion  model  with  a  diffusion  constant 


\  »'■ 


ry******: . 


D 


V2/X. 


In  this  case,  the  two  models  are  said  to  be 


"equivalent" . 

Also,  it  is  argued  in  [Ref.  5]  that  the  detection 
probability  predicted  by  a  constrained  (by  reflecting 
boundaries)  diffusion  model  represents  an  upper  bound  to 
that  predicted  by  the  equivalent  constrained  random  tour 
model.  In  other  words,  the  non-detection  probability 
predicted  by  a  constrained  diffusion  model  is  a  lower  bound 
to  that  predicted  by  the  equivalent  constrained  random  tour 
model.  This  is  reasonable  since  it  is  known  [Ref.  4]  that 
the  target  in  an  unconstrained  diffusion  model  will  "on 
average"  move  a  greater  distance  from  the  origin  than  a 
target  conducting  the  equivalent  random  tour.  Consequently 
the  diffusing  target  would  be  expected  to  encounter  a 
stationary  searcher  more  quickly. 

These  observations,  as  regarding  the  relationship 
between  random  tour  and  its  equivalent  diffusion,  were 
supported  by  plotting  the  results  of  two  simulation 
programs:  RATSIM  and  DIFSIM. 

DIFSIM  (diffusion  simulation)  is  a  Monte-Carlo  search 
simulation  for  a  diffusing  particle  developed  by  Sislioglu 
[Ref.  5]. 

To  generate  the  results  displayed  in  Figure  2.1,  the 

2 

parameters  A  and  R  were  held  fixed  at  1000  nm  and  28.2  nm 


respectively  for  both  programs.  The  diffusion  constant  D 


0000  nm  X(hr”  )  V(nm/hr)  fR?V( 


Figure  2.1  Asymptotic  Approach  of  Random  Tour  to  Diffusion 


used  with  DIFSIM  was  100  nm  /hr.  In  RATSIM  five  different 
( x, V)  pairs  were  selected  such  that  V  /x  =  100  nm  /hr  for 
each  pair.  The  values  of  (x  hr  1/  V  nm/hr)  were  (0.25,  5), 
(1,  10),  (4,  20),  (9,  30),  and  (16,  40). 

It  is  clear  from  Figure  2.1  that,  in  this  case  as  the 
ratio  of  the  characteristic  length  of  the  search  area  to  the 
mean  segment  length  of  the  random  tour  (  fA/V(l/x)  gets 
larger  the  non-detection  probability  curves  for  a  random 
tour  model  asymptotically  approach  that  of  the  equivalent 
diffusion  model. 


C.  MATHEMATICAL  MODEL  OF  DIFFUSION 

Sislioglu  [Ref.  5]  established  a  mathematical  model  to 
predict  the  probability  of  detection  of  a  target  moving 


according  to  the  diffusion  model  (described  in  section  A). 
This  mathematical  model  is  given  by: 


PD(t)  -  1  -  (1  -  ■*!-)  exp  [-  -^5-Pt1 


It  was  later  modified  by  Eagle  [Ref.  3]  to  the  following 
form: 


PD  ( t)  -  1  -  (1  -  •!£-)  exp  [ 


-24.7  RDt 


(A  -  wR  ) 


2.1.5 


where : 


as 


V3:v-» 


m 


Wm 


PD (t)  =  probability  of  detection  at  time  t  in  hr. 

2 

D  =  diffusion  constant,  nm  /hr 

R  =  radius  of  searcher  disk,  nm 

2 

A  =  area  of  search  region,  nm 
So  PND (t)  can  be  given  by: 


PND  ( t)  =  1  -  PD  ( t) 


irR 

PND(t)  =  (1  -  —~)  exp  [ 


-24.7  RDt 


(A  - 


irR2)  1,5 


(2.1 


As  stated  before,  PND(t)  as  given  by  (2.1)  should  represent 
a  lower  bound  on  PND ( t)  as  predicted  by  the  equivalent 
random  tour  model,  and  will  be  used  later  in  the  next 
chapter  as  a  basis  to  derive  the  mathematical  model  of 
random  tour. 

For  simplicity,  equation  (2.1)  will  be  written  in  the 
form: 


PND (t)  *  a  •  e"Bt 


(2.2 


where 


As  indicated  in  [Ref.  b]  ,  the  diffusion  const-u  *-  D  can 

2 

be  approximated  by  V  /x  to  get  a  diffusion  model  equivalent 
to  the  random  tour  with  V  and  x.  So,  if  we  replace  D  by 
V  /X  in  (2.4)  we  get  the  approximate  rate  of  detection  for 
the  equivalent  diffusion  model  in  the  form 


'  r*  i  ;  -  ■ »  ■  m  i.i 


III.  THE  ANALYTICAL  MODEL  FOR  THE  PROBABILITY 
OF  NON-DETECTION 

In  this  chapter  an  experimental  analytical  model  is 
constructed  to  predict  the  probability  of  non-detection  by 
time  t  of  a  target  moving  according  to  the  random  tour  model 
described  in  Chapter  I. 

Simulation  results  from  RAfSlM  will  be  used  as  well  as 
the  relationship  between  the  random  tour  model  and  the 
asymptotically  equivalent  diffusion. 

A.  MODEL  ASSUMPTIONS 

The  following  assumptions  are  made: 

1)  The  target  starting  oosition  is  uniformly  distributed 
over  the  square  searcn  area  A. 

2)  The  target  reflects  perfectly  off  the  area  boundaries. 

3)  The  target  moves  over  the  area  A  according  to  a 
random  tour  with  constant  speed  V  aod  rate  of  course 
change  x. 

4)  The  searcher  is  fixed  at  the  center  of  A. 

5)  The  searcher  detects  with  probability  1  all  targets 
with  a  range  of  R  or  less.  The  searcher  never  detects 
targets  at  ranges  greater  than  R.  (That  is,  the 
searcher  has  a  "cookie-cutter"  sensor  with  detection 
range  R  [Ref.  2 ] . ) 

6)  The  problem  ends  when  the  target  is  detected. 

B.  CLASSIFICATION  OF  VARIABLES 


1.  The  Independent  Variables 

2 

Search  area  A  in  square  nautical  miles  (nm  ) 


•  Probability  of  non-detection  by  time  t,  PND(t), 
i.e. ,  PND ( t)  =  f (A,  V,  R,  X)  . 

C.  CONSTRUCTION  OF  THE  MODEL 

By  plotting  PND(t)  versus  t,  as  estimated  by  RATSIM  and 
with  a  logarithmic  scale  for  the  Y-axis,  it  was  observed 
that  the  resulting  curves  were  very  nearly  linear  with 
negative  slopes  (see  Figure  3.1). 

This  linear  relationship  on  a  logarithmic  scale  graph 
suggests  the  following  functional  form  for  PND(t): 

PND (t)  =  a  •  e"Yt  .  (3.1 

In  the  course  of  this  research  approximately  300 
simulation  experiments  with  RATSIM  were  conducted.  All 
showed  PND ( t)  to  be  approximately  given  by  (3.1).  Figures 
3.2  through  3.5  are  representative. 

This  thesis  attempts  to  fit  the  simulation  data  and 
establish  values  of  a  and  -y  as  functions  of  the  problem 
independent  variables  A,  V,  R  and  X. 

A  small  subroutine  was  added  to  the  main  program  of 
RATSIM  to  compute  a  least-squares  estimate  of  a  and  y.  The 
formulas  used  [Ref.  6]  were: 


i  t  HOUUCED  AT  Gu  .  .  itl.IFNI  FXPF.NOE 


^vtv,v,v;v.v,Y, 


— r*.  - 

*>w 


fc.PRODLK.ED 


As  would  be  expected,  all  simulations  conducted  showed 


PND (0  +  )  -  1  -  wR2/A 


2.  Submodel  for  y 

As  stated  before,  PND (t)  predicted  by  a  diffusion 
model  appears  to  be  a  lower  bound  to  those  values  predicted 
by  the  equivalent  random  tour  model. 

A  study  of  the  simulation  data  suggests  that  y  in 
equation  (3.1)  can  be  estimated  by: 


y  »  6(1  -  e"*) 


(3.3) 


where 


is  the  detection  rate  of  the  equivalent 
diffusion  model  given  by  equation  (2.5). 

is  a  function  of  the  independent  problem 
variables  A,  V,  R  and  x. 


In  this  thesis  an  attempt  is  made  to  find  the 
functional  relationship  between  *  and  the  independent 
variables.  It  should  be  noted  that  there  may  be  other 
functional  forms  for  y  that  fit  the  simulation  data  as  well 
or  better  than  equation  (3.3). 


»< 


3.  Submodel  for  » 

This  submodel  includes  the  A,  R,  V  and  x.  So  it  can 
be  expected  to  be  more  complex  than  the  submodel  for  ct.  To 
simplify  the  problem,  the  relationship  between  <i>  and  each 
one  of  these  variables  was  investigated  separately  at  first. 
Then  a  combination  of  these  separate  relationships  was  used 
to  construct  the  required  submodel  for  4>. 

a.  Relationship  Between  *  and  Area  Size  A 

To  obtain  this  relationship,  the  variables  V,  R 
and  x  were  held  fixed  at  10  nm/hr,  10  nm  and  1  hr  ^ 

2 

respectively.  The  area  size  A  was  varied  between  900  nm 

2 

and  20000  nm  . 

Each  simulation  run  required  the  independent 
variables  A,  R,  V  and  x  to  be  specified,  and  gave  a  best  fit 
for  t  as  an  output.  Then  from  equation  (3.3)  we  have 

♦  -  -  In  (1  -  *  )  •  (3.4) 

Substituting  (2.5)  into  (3.4)  yields 

+  «  -In  {1  -  T*<a  *r  > -  )  .  (3.5) 

24.7  Rv 

By  plotting  the  values  of  tr  calculated  by 


equation  (3.5),  versus  the  corresponding  values  of  A,  it  was 
found  that  a  power  function  fit  the  data  very  well  (see 


Figure  3.6).  The  least-squares  best-fit  4>  was  found  to  be 
given  by 


*  =  0. 0080678 (A) °*  5028 


This  implies  that 


*  « 


(3.6) 


where 


"«"  means  "is  proportional  to". 

b.  Relationship  Between  *  and  Target  Speed  V 

Here  the  variables  A,  R  and  *  were  held  fixed  at 
2  -1 

10000  nm  ,  10  nm,  and  1  hr  respectively.  V  was  then 
varied  between  2.5  nm/hr  and  25  nm/hr.  The  simulation 
output  y  and  equation  (3.5)  were  used  to  generate  the 
corresponding  values  of  *. 

By  plotting  *  versus  V  and  fitting  a  power 
function  to  the  data,  it  was  observed  that 

*  -  8. 3357  (V)-1,001 


(See  Figure  3.7).  This  implies  that 


HNMFN1  f.  .xrf  NSE 


o 


PRODUCED  AT  Go.  ..i.IENI  EXPENSE 


c.  Relationship  Between  *  and  the  Rate  of  Course 
Change  x 

Now  A,  V  and  R  were  held  fixed  at  10000  nm2, 

10  nm/hr  and  10  nm  respectively.  Then  x  assumed  the 
following  values:  0.01,  0.05,  0.1,  0.2,  0.5,  1,  1.5,  2, 
2.5,  3,  and  4  hr-1.  The  resulting  best  fit  power  function 
(see  Figure  3.8)  was  found  to  be 

+  =  0. 83419  (X)  1*0003 

This  means  that 

*  -  X  (.38 

d.  Relationship  Between  *  and  Detection  Radius  R 

The  following  values  were  assigned  to  R:  2.5, 

5,  7.5,  10,  12.5,  15,  20,  25,  30,  and  35  nm.  The  other 

o 

variables  A,  V  and  x  were  held  fixed  at  10000  nm  ,  10  nm/hr 
and  1  hr-1  respectively.  The  best  fit  power  function  was 

*  =  0.8423(R)-0,009 

This  indicates  that  *  is  nearly  independent  of  R  over  this 
range  of  R  values  (see  Figure  3.9). 


e.  Summary  and  Conclusions 


Now,  we  can  summarize  the  previous  relationships 


as  follows: 


4»  «  A  f 


*  «  1/V 


iji  «  X 


This  leads  to  the  following  conclusion: 


ijj  oc 


T7T  x 


Or  equivalently, 


V 


(3.9) 


where  K  is  a  proportionality  constant  to  be  estimated 
from  the  simulation  data. 

f.  Estimation  of  the  Coefficient  K 

The  outputs  y  of  156  simulation  experiments  with 
RATSIM  were  used  to  produce  146  sample  K  values.  The  value 
of  K  was  calculated  from  equation  (3.9)  as  follows: 


fc-JW 

s-V-'Cv 

v'.\W 
c  *.W 


!>\*  V 


A  *\  •• 

>.v:> 


i>v 

fA  X 


where  the  value  of  *  was  determined  by  equation  (3.5) 


(3.10) 


The  156  RATSIM  experiments  used  to  estimate  K 
resulted  in  a  sample  mean  o£  0.084  and  a  sample  standard 
deviation  of  0.0016.  These  data  suggest  that  with 
probability  0.9,  K  lies  in  the  interval  [0.082,  0.087].  Th 
bounds  of  this  confidence  interval  were  the  observed  5  and 
95  percentile  points. 

The  histogram  and  the  statistical  summary  table 
for  this  data  are  displayed  in  Figure  3.10. 
g.  Final  Submodel  for  * 

By  estimating  the  value  of  K  and  applying 
equation  (3.9)  we  can  construct  the  final  submodel  for  *  as 
follows : 


*  *  0.084 


nr  x 
1  ..  ■  ■  ♦ 


(3.11 


h.  Final  Submodel  for  y 

Substituting  equations  (3.11)  and  (2.5)  into 


equation  (3.3)  we  get 


24.7  RV  r.  ,  _  A_.  fX 

o  i  ?  [1  “  exp(-0.084  w  )]  • 

(A  -  irR  )  *5X  V 


(3.12 


4.  The  Final  Form  of  the  Random  Tour  Analytical  Model 
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IV.  VERIFICATION  OF  THE  RANDOM  TOUR  MODEL 


I 


A.  DIMENSIONAL  ANALYSIS 

From  equation  (3.3)  it  is  clear  that  *  must  be 
dimensionless.  Now  writing  l  as  K  fA~  x/V,  we  see  that  K  also 
is  a  dimensionless  coefficient 

B.  LIMIT  OF  y  AS  X  ♦  0 

From  equation  (3.12)  we  have 


Lim  y  -  Lim  {  [1  -  exp(-0.084  &^)  ]  } 

X-0  X-0  (A  -  wR  ; a 


24.7  RV2 
(A  -  irR2)1, 


5 


(0.084 


2.0748  RV 
(A  -  wR2) 


SK 


1.5 


(4.1) 


So  if  equation  (3.12)  is  a  reasonable  estimate  for  y, 
then  as  x  ♦  0  RATSIM  should  give  a  best  fit  y  given  by 
equation  (4.1).  To  test  this,  four  groups  of  simulation 
experiments  were  conducted.  In  each  group,  values  of  A,  R 
and  V  were  held  constant  and  x  was  varied  from  10  to  0.01. 
Figure  4.1  shows  the  best  fit  y  plotted  against  1/x  for 
each  of  the  simulation  groups.  Also  plotted  is  a  horizontal 


line  intersecting  the  Y-axis  at  the  value  given  by  equation 
(4.1). 


For  these  simulations,  it  appears  that  equation  (3.12) 
holds  as  x  *  0. 

It  is  noted  that  X  =  0  means  that  the  target  never 
changes  course  except  when  reflecting  off  the  area 
boundaries. 

C.  LIMIT  OF  y  AS  *R2  -  A 

2 

From  equation  (3.12)  we  see  that  y  ♦  •  as  *R  ♦A, 
which  implies  that  PND(t)  ♦  0  for  t  >  0.  This  is  as  would 
be  expected. 

D.  ASYMPTOTIC  APPROACH  TO  DIFFUSION  MODEL 

As  stated  before,  when  the  ratio  of  the  characteristic 

length  of  the  search  area  to  the  mean  segment  length  of  the 

random  tour  (  )|~a7v(1/X))  becomes  large,  we  find  that  the 

random  tour  model  approaches  the  asymptotically  equivalent 

2 

diffusion  model  with  a  diffusion  constant  V  /X.  This  is 
consistent  with  equation  (3.12).  Since  by  taking  the  limits 
of  both  sides  of  (3.12)  as  (  )[a/V(1/x))  -  we  get 


Lim  y 
IJEa/V-*— 


24.7  RV2 
(A  -  wR2)1,5X 


This  is  the  value  of  0  given  by  equation  (2.5)  for  a 

2 

diffusion  model  when  the  diffusion  constant  is  V  /x. 


E.  ASYMPTOTIC  APPROACH  TO  RANDOM  SEARCH  MODEL 


The  random  tour  model  of  Koopman  [Ref.  7]  predicts  that 
the  detection  rate  of  a  randomly  moving  target  is  2RV/A. 

As  fX  a/v  *  0,  the  model  presented  here  results  in  a 
detection  rate  of 


2.0748  RV  VA 
(A  -  .R2)3/2 


2 

For  small  »R  /A  these  two  expressions  are  nearly  equal. 

F.  LIMIT  OF  y  AS  V  ♦  0 

By  taking  the  limits  of  both  sides  of  equation  (3.12)  as 
V  ♦  0  we  get 


Lim  Y  =  0 
V*0 


This  means  that  for  V  -  0, 

PND(t)  »  1  -  irR2/A  ,  t  >  0  , 


which  is  as  would  be  expected. 

G.  SENSITIVITY  ANALYSIS 

Figure  4.2  illustrates  how  equation  (3.13)  behaves  as 
the  independent  variables  A,  V,  R  and  A  are  varied  one  at  a 
time.  The  base  case  considered  was: 


A  =  10,000  nm  , 


V  *  10  nm/hr, 

x  *  1  hr 

R  =  10  nm, 

t  =  20  hr 

Equation  (3.  13)  is  seen  to  be  an  increasing  function  of 
A  and  x,  and  a  decreasing  function  of  V  and  R.  This  agrees 
with  intuition. 

As  A  increases,  the  target  has  more  area  in  which  to 
hide.  So  PND  will  increase. 

As  V  or  R  increases,  the  target  will  be  more  likely  to 
encounter  the  detection  disk.  So  PND  decreases. 

And  as  x  increases,  the  target  tends  to  remain  closer  to 
its  starting  position.  So  PND  will  increase. 

H.  FINAL  VERIFICATION 

There  exist  no  actual  data  available  from  real  life 
observations.  Therefore,  the  output  of  RATSIM  was  used  for 
final  verification  of  the  model. 

To  achieve  this  purpose  47  combinations  of  different 
values  of  the  independent  variables  A,  V,  R  and  X  were  used 
as  input  to  both  simulation  program  RATSIM  and  the  proposed 
analytical  model  given  by  equation  (3.13). 

These  47  experiments  were  classified  into  four  groups, 
where  in  each  group  only  one  parameter  was  varied  while  the 

50 


2 

others  were  kept  at  the  base  case  value  (A  =  10000  nm  ,  V  = 
10  nm/hr,  R  *  10  nm  and  X  =  1  hr  *) .  The  outputs  of  these 
different  experiments  are  displayed  in  Table  4.1. 

By  looking  carefully  into  the  values  displayed  in  Table 
4.1,  we  observe  that  there  is  a  little  difference  between 
the  values  obtained  from  simulation  and  the  corresponding 
values  estimated  by  the  proposed  analytical  model,  except 
for  large  values  of  x  (X  >  20),  and  for  large  values  of 
wR2/A  ( wR 2/A  >  0.3). 

So,  we  can  say  that  the  proposed  analytical  model  is 
reasonable  for  the  realistic  values  of  the  problem 
independent  variables  (A,  V,  R  and  x)  used  in  antisubmarine 
warfare  (ASW) . 

Figures  4.3  and  4.4  show  a  comparison  of  PND ( t) 
generated  by  RATSIM  and  the  analytical  model  for 
representative  values  of  the  independent  variables.  For 
many  cases  the  fit  is  so  close  that  the  curves  are  nearly 
Indistinguishable. 


TABLE  4.1 

VARIATION  OF  y  AND  (MAX  A|*  WITH  VARIATION  IN  THE 
INDEPENDENT  VARIABLES 


(1) 
Vary  A 


(2) 

Model  y 


(3) 

Simulated  y 


(4) 

| Max  A | 


=  400  nm 

900 
2000 
4000 
6000 
8000 
10000 
12000 
14000 
16000 
18000 
20000 
25000 


Vary  R 


4.8 
0.  388 
0.112 
0. 0455 
0.0276 
0. 0194 
0.0147 
0. 0118 
0.00972 
0. 00823 
0.0071 
0. 00622 
0.00468 


2.587 
0.418 
0.117 
0. 0461 
0.0278 
0. 0193 
0.0145 
0. 0117 
0.00975 
0. 00825 
0.00721 
0. 00625 
0.00466 


0.08 
0.0305 
0.0301 
0.0225 
0.0165 
0.  011 
0.012 
0.01 
0.015 
0.013 
0.0123 
0.015 
0.016 


.  5  nm 

0.00352 

0. 00398 

0.  047 

5 

0.0072 

0.00738 

0.0345 

10 

0. 0147 

0. 0145 

0.012 

15 

0.0235 

0.0229 

0.013 

20 

0. 0343 

0.0338 

0.011 

25 

0.0487 

0.048 

0.017 

30 

0. 0693 

0. 0683 

0. 0216 

35 

0.1 

0.0986 

0.032 

40 

0.  16 

0. 15 

0.  041 

50 

0.706 

0.308 

0.07 

|Max  a|:  The  maximum  absolute  difference  between  PND(t) 
estimated  by  simulation  and  PND(t)  estimated  by 
the  analytical  model  at  the  same  t,  over  the 
whole  experiment  period  (TMAX) . 
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TABLE  4.1  (CONTINUED) 


(2) 


(3) 


(4) 


0.00156 

0. 00177 

0. 0323 

0.00527 

0.00526 

0.0244 

0. 0147 

0.0145 

0.  012 

0.025 

0.0251 

0.014 

0. 0355 

0. 0356 

0.  012 

0.0462 

0.0472 

0.013 

0.  057 

0. 0582 

0.  021 

0.06773 

0.0689 

0.014 

0. 0785 

0. 0801 

0.  019 

0.1 

0.11 

0.0305 

0. 2088 

0. 2205 

0. 0329 

0.4268 

0.446 

0.017 

0. 0209 

0. 0224 

0.0217 

0.02004 

0.0217 

0.0162 

0. 0184 

0. 0196 

0. 0154 

0.017 

0.0168 

0.014 

0. 0158 

0.01579 

0. 0092 

0.0147 

0.0145 

0.012 

0. 0137 

0. 0125 

0.  016 

0.0105 

0.0102 

0.0242 

0.0051 

0. 00512 

0.  037 

0.00259 

0.00294 

0.0521 

0.00173 

0. 00224 

0.  054 

0.0013 

0.00687 

0.065 
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PROBABILISTIC  ANALYSIS  OF  THE  MODEL 


A.  CUMULATIVE  DISTRIBUTION  FUNCTION  (CDF) 

Let  T  be  the  random  variable  for  time  of  detection.  An 
let  F(t)  be  the  cumulative  distribution  function  (CDF)  for 
T.  That  is, 

P  {T  <  t}  =  F  (t)  • 

The  model  presented  here  implies  that  F(t)  can  be 
closely  approximated  by 

F (t)  =  u(t)  (1  -  ae~ytJ  •  (5.1 


where 

u(t)  is  0  for  t  <  0  and  1  for  t  >  0 

It  is  noted  that  equation  (5.1)  satisfies  the  following 
properties  of  a  CDF: 

1)  Lim  F(t)  *  1, 

2)  F  (0)  -  0, 

3)  F(t)  >  0, 

4)  F ( t)  is  a  non-decreasing  function  (see  Figure  5.1). 


B.  DENSITY  FUNCTION 

By  taking  the  first  derivative  of  F(t)  with  respect  to 
time,  we  can  derive  the  density  function  (f(t))  for  T  as 
follows : 


f(t)  ■  dF(t)/dt 

-  u { t )  (aye'Yt)  +  6 (t) (1  -  a)  . 


(5. 


where 


«(t)  is  the  Dirac  «-function. 

C.  EXPECTED  VALUE  OF  DETECTION  TIME 

The  expected  detection  time  E[T]  can  be  derived  as 
follows : 


E[T]  *  J  [i  -  F  (t)  ]dt 
0 

m 

»  /  [1  -  u(t)  (1  -  ae'Tt))dt 


F 


Replacing  a,  y  by  their  expressions  given  by  equations 
),  (3.12)  respectively  we  get: 


24.7  RV  [  (exp  (-0. 084  -^)  )  -  1] 

(5.4 

This  equation  shows  how  E[T]  varies  with  the  problem 
independent  variables  A,  R,  V  and  X.  The  variation  of  E[TJ 
with  each  of  these  variables  is  indicated  in  Figure  5.2. 


D.  CONDITIONAL  CDF 

If  we  assume  that  there  will  be  no  detection  at  the 
beginning  of  the  search  period,  we  may  derive  the  following 
conditional  CDF  (FQ(t))s 


F0(t)  =  P  {Detection  by  time  t  I  no  det.  at  t  =  0+) 


*  P  {T  <  t  |  T  >  0} 


P  (T  >  0,  T  <  t] 
P  (T  >  0) 


P  (0  <  T  <  t] 
P  (T  >  0) 


If  we  substitute  t  *  0+  in  (5.1)  we  get 
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P  (0  <  T  <  t)  =  F  (t)  -  F  (0  ) 


(1  -  ae"Yt)  -  (1  -  a) 


=  a(l  -  e_Yt)  . 


Also , 


P (T  >  0)  =  F(0+) 


=  1  -  F  (0  ) 


Substituting  (5.7)  and  (5.8)  into  (5.5),  we  get 


FQ(t)  *  (1  -  e"Yt)  . 


This  function  (5.9)  is  a  CDF  for  an  exponential 
distribution  with  parameter  y  (an  expression  for  y  is  given 
by  (3.12)). 


E.  CONDITIONAL  EXPECTED  VALUE  OF  DETECTION  TIME 

The  conditional  expected  first  detection  time  E[Tq]  can 


be  defined  as  follows: 


E[Tn]  ■  E[T  |  no  detection  at  t  >  0] 


-  f  F0{t)  dt 
0 


1 

~  a 

Y 


(5.10) 


F.  CONDITIONAL  DENSITY  FUNCTION 

Finally,  the  conditional  density  function  ( f 0 ( t ) )  can  be 
derived  as  follows: 


fott) 


dF0(t) 

3t 


Y 


e-Yt 


(5.11) 


If  we  compare  (5.3)  and  (5.10),  we  will  observe  that 

—  <  —  ,  since  o  <  1  for  R  >  0. 

Y  Y 

This  implies  that  the  conditional  expected  first  detection 
time  is  greater  than  the  unconditional  one.  This  is 
reasonable,  since  in  the  unconditional  case  we  have  an 
opportunity  to  detect  the  target  at  time  0  .  This 
conclusion  is  demonstrated  clearly  by  comparing  F(t)  and 
Fq { t)  as  illustrated  in  Figure  5.3,  where  we  always  find 
that  Fq (t )  is  less  than  F(t)  at  any  value  of  t,  except  at 
t  ■  •  where  F(t)  ■  FQ(t)  -  1. 
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APPENDIX 


RATSIM  COMPUTER  PROGRAM 

In  order  to  give  access  to  the  logic  used  in  building 
the  simulation  model  RATSIM,  a  complete  program  listing  is 
included  in  this  appendix  following  the  list  of  variables 
used  in  the  simulation  model. 

LIST  OF  VARIABLES 

The  variables  used  in  the  simulation  model  are  listed 

below  according  to  their  first  appearance  in  the  program: 

R  =  Radius  of  searcher  detection  disk  in  nautical 

miles 

V  ■  Speed  of  target  in  nautical  miles  per  hour 

Inc  *  Time  increment  for  each  discrete  step  in 

minutes 

REP  *  Number  of  replications 

TMAX  *  Detection  period  in  minutes 

SUM (I)  -  Number  of  detections  at  time  increment  I 

HIST (I)  *  Accumulative  number  of  detections  up  to 
increment  I. 

POSX  *  X  component  of  target's  position 

POSY  *  Y  component  of  target's  position 

XS  ■  X  component  of  target's  starting  position 

YS  ■  Y  component  of  target's  starting  position 

*  Course  9  in  radians 


ANG 


TLEG 


Time  leg  for  each  segment 


D  ®  Distance  between  the  target  location  and  the 

center  of  the  detection  disk 

XN  =  X  component  of  the  target's  new  position 

YN  *  Y  component  of  the  target's  new  position 

PROBD  =  Probability  of  detection 

PROBS  =  Probability  of  non-detection 

A  =  1  -  itR2/A  (a) 

B  =  Detection  rate  (y) 
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*★★★***★******★*★*★★★*★**★★★******★★*★**★****★**★******** 
* PROGRAM  NAME:  RATSIM  * 

♦THIS  PROGRAM  SIMULATES  2-DIMENSIONAL  RANDOM-TOUR  MODEL  * 

REAL  INC, L,  L2 , LAMBDA 

INTEGER  REP ,CTR, TIME , SUM( 8000) ,  HIST(8000) 

DIMENSION  XS ( 2500 ), YS ( 2500 ) , TEXP ( 3000 ) , TH( 3000) 

DIMENSION  TI ( 8000 )  , PROBD( 8000) , PROBS ( 8000) 

DIMENSION  Y ( 8000 )  ,  Z  ( 8000 ) 

DOUBLE  PRECISION  DSEED 
DSEED=89456 . DO 
NR=2400 

CALL  GGUBS  (DSEED, NR,  XS) 

DSEED=73452 . DO 
NR=2400 

CALL  GGUBS  (DSEED, NR,  YS) 

R=10. 

AREA= 10000. 

V=10. 

LAMBDA=1 . 

INC=3 . 

REP=  2400 
TMAX=100*60 
L=AREA*  * . 5 
L2=L*2 
SER=L/2 . 

MAXCTR= I NT ( TMAX/ I NC ) + 1 
DO  10  1=1, MAXCTR 

SUM( I )=0 
HIST( I )=0 
CONTINUE 
DO  50  1=1, REP 

DSEED=6095 . DO*DBLE ( FLOAT ( I ) ) 

NR=3000 

CALL  GGUBS  (DSEED, NR, TH) 

DSEED=2211 . 0D0*DBLE( FLOAT ( I ) ) 

XM= 1/LAMBDA 
NR=3000 

CALL  GGEXN (DSEED , XM , NR , TEXP ) 

POSX=XS( I  )*L 
POSY=YS( I  )*L 
TIME=0 
CTR=1 

DO  40  J=1 , 3000 

ANG=6 . 2832  *TH(J) 

TLEG=  TEXP (J) *60 
N=INT(TLEG/INC) 

DO  35  M=1,N 

D=( ( (POSX-SER)**2)+( (POSY-SER) **2 ) ) ** . 5 

CHECK  FOR  DETECTION 

IF(D.LE.R)  GO  TO  45 

IF  (CTR.GT. MAXCTR)  GO  TO  48 

CURRENT  POSITION  OF  THE  TARGET 

XN=POSX+V* ( INC/60 ) *S IN ( ANG ) 

YN=P0SY+V* ( INC/60 ) *COS ( ANG) 

66 


m 

Sf&s 
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CHECK  FOR  REFLECTION 
IF(XN)  16,16,17 
XN=-XN 

ANG=6 . 2832-ANG 

IF(XN-L)  19,18,18 

XN=L2-XN 

ANG=6 . 2832  -ANG 

IF(YN)  20,20,21 

YN=-YN 

ANG=3 . 1416-ANG 
IF(YN-L)  30,22,22 
YN=L2-YN 
ANG=3. 1416-ANG 
POSX=XN 
POSY=YN 
CTR=CTR+1 
TIME=TIME+ INC 
CONTINUE 
CONTINUE 

INDEX=INT( TIME/INC) +1 
SUM ( INDEX )=SUM( INDEX) +1 
CONTINUE 
CONTINUE 
HIST( 1 )=SUM( 1 ) 

PROBD ( 1 ) =FLOAT ( HI ST ( 1 ) ) /FLOAT ( REP ) 

PROBS (1)=1  -PROBD ( 1 ) 

TI ( 1 )=0 . 0 

DO  300  I =2 , MAXCTR 

HIST( I )=SUM( I ) +HIST( 1-1) 

PROBD ( I ) =FLOAT ( HI ST ( I ) ) /FLOAT ( REP ) 

PROBS( I )=1  -PROBD ( I ) 

TI( I )=FLOAT(I-l)*( INC/60) 

WRITE( 6, 250)  TI ( I ) , PROBS ( I ) 

FORMAT(2X, F6. 2, 2X, F16.ll) 

CONTINUE 

************************************************* 
THE  FOLLOWING  SUBROUTINE  IS  TO  ESTIMATE  THE  DETECTION 
RATE"B"  AND  THE  COEFFICIENT  "A"  BY  USING  A  NUMERICAL 
METHOD . 

DO  350  1=1, MAXCTR 
Y(I)=0 
2(1) =0 
CONTINUE 

DO  380  1=1, MAXCTR 
IF( PROBS ( I ) . LE . 0 )  GO  TO  390 
Y( I )=PROBS( I ) 

Z( I )=TI ( I ) 

CONTINUE 
K=I 
SZ=0 
SZ2=0 
SLY=0 
SLYZ=0 
SLY2=0 

DO  400  J=1,K 


v 


SZ=SZ+Z( J) 

SZ2=SZ2+Z( J)**2 
IF(Y(J) .LE.O)  GO  TO  500 
SLY=SLY+ ALOG ( Y ( J ) ) 
SLYZ=SLYZ+Z ( J ) * ALOG ( Y ( J ) ) 
SLY2=SLY2+ALOG(Y(J) )**2 
CONTINUE 

Ul=( SZ2*SLY) - ( SZ*SLYZ ) 
G=(K*SZ2 ) - ( SZ**2 ) 

U2= ( K*  SLYZ ) - ( SZ*  SLY ) 

F=U1/G 

B=U2/G 

A=EXP(F) 

WRITE ( 6, * )  A, B 


LIST  OF  REFERENCES 


Washburn,  A.  R. ,  Search  and  Detection,  Military 
Applications  Section,  Operations  Research  Society  of 
America,  1981. 

Washburn,  A.  R. ,  "Probability  Density  of  a  Moving 
Particle,"  Operations  Research,  V.  17,  1,  September 
1969. 


Eagle,  J.  N. ,  Report,  Estimating  The  Probability  of 
Diffusing  Target  Encountering  a  Stationary  Sensor.  J 


985. 


etection  Probability  of  a  Diffusion  Target,  Master's 
Thesis,  Naval  Postgraduate  School,  Monterey,  California 
September  1984. 

Dorn,  William  S.  and  McCracken,  Daniel  D. ,  Numerical 
Methods  with  Fortran  IV  Case  Studies.  John  Wiley  &  Sons 


Office  of  the  Chief  of  Naval  Operations,  Navy 
Department,  Washington,  D.C.,  Operations  Evaluation 
Group  Report  No.  56,  Search  and  Screening,  by  B.  0. 
Koopman,  1946. 


j :  ,V-V.V-  'V^V  -  ^  ~  " 


INITIAL  DISTRIBUTION  LIST 


Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  Virginia  223Q4-6145 

Library,  Code  0142 

Naval  Postgraduate  School 

Monterey,  California  93943-5100 

Professor  James  N.  Eagle,  Code  55Er 
Naval  Postgraduate  School 
Monterey,  California  93943-5100 

Professor  A.  R.  Washburn,  Code  55 
Naval  Postgraduate  School 
Monterey,  California  93943-5100 

Professor  Edward  B.  Rockower,  Code  55Rf 
Naval  Postgraduate  School 
Monterey,  California  93943-5100 

Abdel-Aziz  Al-Bassiouni 
SMC  Box  1689 

Naval  Postgraduate  School 
Monterey,  California  93943-5100 

COL  Aref,  Salama 
Operations  Research  Center 
Ministry  of  Defense 
Cairo,  Egypt 

Salah  Ibrahim  Abd  El-Fadeel 

16,  Samaan  St  -  Guizirat  Badran  St., 

Shoubra 

Cairo,  Egypt 


2 


